Melting of generalized Wigner crystals in transition metal dichalcogenide heterobilayer Moiré systems

Moiré superlattice systems such as transition metal dichalcogenide heterobilayers have garnered significant recent interest due to their promising utility as tunable solid state simulators. Recent experiments on a WSe2/WS2 heterobilayer detected incompressible charge ordered states that one can view as generalized Wigner crystals. The tunability of the transition metal dichalcogenide heterobilayer Moiré system presents an opportunity to study the rich set of possible phases upon melting these charge-ordered states. Here we use Monte Carlo simulations to study these intermediate phases in between incompressible charge-ordered states in the strong coupling limit. We find two distinct stripe solid states to be each preceded by distinct types of nematic states. In particular, we discover microscopic mechanisms that stabilize each of the nematic states, whose order parameter transforms as the two-dimensional E representation of the Moiré lattice point group. Our results provide a testable experimental prediction of where both types of nematic occur, and elucidate the microscopic mechanism driving their formation.

We are interested in the case of a fixed number of particles M ≤ N , i.e. the only valid particle configurations Λ are those with cardinality M . We treat elements of the point group τ ∈ G as maps τ : Z N → Z N which map lattice site indices to lattice site indices. Particle exchange is a map η : Z N × Z N × G → Z N defined as Note that η preserves the cardinality of Λ.
The Hamiltonian for the system is H(Λ) = 1 Algorithm: 1. Fix a particle configuration Λ, and initialize an empty set C (the cluster) 2. Choose an order-2 element, τ * , of the lattice point group.
3. Randomly choose a site i, set C = C ∪ {i, τ * (i)} and Λ = η(Λ, i, τ * ) 4. For each other site k with V ik = 0 and k / ∈ C: (a) Calculate The form of ∆ ik is chosen to ensure detailed balance (more detail in the next section).
(b) With probability max 0, 1 − e −β∆ ik (Λ) , set C = C ∪ {k, τ * (k)}, Λ = η(Λ , k, τ * ), and record k in a stack data structure 5. Pop an element j from the stack and repeat step 4 with j playing the role of i 6. Repeat step 5 until the stack is empty and then return the updated particle configuration Λ .

Proof of Detailed Balance:
To prove that this generates the correct equilibrium probability distribution (Boltzmann distribution), we need to show that the probability of a particle configuration Λ moving to a configuration Λ , P(Λ → Λ ), satisfies the detailed balance condition, i.e.
Observe that, for a move corresponding to a cluster C, we can write P(Λ → Λ ) = P in (C, Λ)P out (C, Λ) where the first factor is the probability of forming the cluster containing the sites in C and the second is the probability that no sites in Z N \ C are included in C. Similarly, we write P(Λ → Λ) = P in (C, Λ )P out (C, Λ ). Because (1) τ * is order-2 and a symmetry of H and (2) P in (C, Λ) only depends on lattice sites in C, we have that P in (C, Λ) = P in (C, Λ ).
By construction of the algorithm, we have (using standard probability rules) that and where in the last equation we have used the fact that for i ∈ C, n(Λ , i) = n(Λ, τ * (i)). Thus we have that where in the second equality we have used that τ * is order-2 and a symmetry of the Hamiltonian. Note also the factor of 1/2, which is due to the fact that since i ∈ C ⇔ τ * (i) ∈ C and ∆ ik (Λ) = ∆ τ * (i)k (Λ) the sum double counts. What remains is to show that the RHS of eq. 1 is the same as eq. 2. Using the fact that τ * is order-2 and a symmetry of the Hamiltonian, we can write It is then easy to see that the difference This completes the proof and establishes that the algorithm defined above satisfies detailed balance.
A Note about Ergodicity: Note that by choosing τ * as an appropriate reflection, it is possible to form a cluster consisting of a nearest neighbor particle exchange with finite probability. Therefore, it is possible with finite probability to get from any particle configuration to any other in the same way that it is possible to realize any spin configuration in an Ising model via a sequence of single spin flips. Thus we can see that this algorithm is indeed ergodic.  Fig. 1. Wire Energies. Energy of isolated, continuous wires interacting via eq. ?? with uniform line charge density λ as a function of wire length L. We show results for parallel wires seperated by 5a/2 (purple) and wires angled at π/3 with minimum separation √ 3a/2 (green). For sufficiently long wires the angled wires have lower energy.

SECTION 2: CALCULATION OF NEMATIC CORRELATION FUNCTION AND ORIENTATION
To calculate the nematic order parameter expectation value N ( r) = n( r)e i2θ( r) in Monte Carlo, we write N ( r) as where the vectors δ point to the nearest neighbors of r and ρ( r) is the occupation number of the site at r. This is easy to calculate directly from a Monte Carlo configuration and average over the Markov chain. To calculate the orientation cos(6θ) we use the fact that we can write the vector part of the nematic order parameter as (cos(2θ), sin(2θ)) T . From this we can simply calculate θ and hence cos(6θ) from the nematic order parameter and average it over the Markov chain.

SECTION 3: ENERGETICS OF π/3 INTERSECTION AND PARALLEL DIMER COLUMN FRAGMENTS
We can understand why the π/3 intersections are energetically preferred schematically by calculating the energy for two continuous, isolated wires of uniform line charge density λ interacting via eq. ?? as a function of wire length L. We approximate the wires as being at the center of the dimer columns, and plot the results for parallel wires separated by 5a/2 in the purrple curve in Fig. S1. For a π/3 intersection with wires at the center of the dimers, the wires terminate with a separation of √ 3a/2. We show the results for such wires in the green curve in Fig. S1. Indeed, the energy of the parallel wires exceeds that of the angled wires for sufficient L.

SECTION 4: THE SYSTEM SIZE DEPENDENCES
Here we compare simulations on the = 20 lattice to simulations on a smaller, = 10 lattice. At = 10, the 2π/3 domain wall junction appears at ν = 0.36 in Fig. 2(a). This is the defining feature of the hexagonal domain wall state. However, the system is still not large enough to fit an entire hexagonal domain. Increasing the density to ν = 0.38 in Fig. 2(b), we can see that π/3 intersections begin to emerge as we would expect for the nematic-II state, but the system is still too small to fit longer columnar fragments. Finally, in Fig 2(c) at ν = 0.48, the simulation cell isn't big enough to fit stripe fragments and dislocations in a single orientation, and instead forms two domains of stripe fragments. What is encouraging about the results from the = 10 system is that the defining features of the  Fig. 2. Results from = 10 system. Low-temperature Monte Carlo snapshots from an = 10 system at densities corresponding to the intermediate phases observed in the larger = 20 system (c.f. main text Fig. 1(d)). (a) ν = 0.36 corresponding to the hexagonal domain wall state (b) ν = 0.38 corresponding to the type-II nematic state (c) ν = 0.48 corresponding to the type-I nematic state.
intermediate states show up as long as they fit in the smaller system. However, if the system size is too small for a particular structure to show up (e.g. dislocations), it may not be visible.